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Problem: 

Given  the  speed-of-sound  distribution  in  the  sea  the  problem  is  to 
calculate  the  trace  of  a sound  lay  which  loaves  a prescribed  location  in  a 
prescribed  direction.  The  formulation  includes  the  curvature-of-the-earth 

1 a 

factor  and  other  curvature  terms  arising'in  the  selection  of  horizontal 
coordinates.  K 

Part  A:  Formulation  of  the  Equations 

].  The  Governing  Vector  Equation 

The  ray  is  directed  normal  to,  and  in  the  direction  of  advancement 
of,  the  wave  front.  We  assign  the  unit  vector  t “ to  this  direction.  At  a 
point  in  the  progression  of  the  sound  ray  the  unit  vector  ^ and  the  vector 
Vc  , being  the  local  ascendent  of  the  sound  speed,  together  generally 
define  a plane.  This  plane  is  the  plane  of  Pig.  1. 


Pig  . 1 Refraction  of  the  Ray 


l0r  arrives  at  (The  integration  may  also  be  performed  backwards  in  time.) 
^Geometric  vectors  are  denoted  by  curl  underscore. 


The  component  of  Vc  which  Is  normal  to  t , and  hence  lies  in  the 
wave  front,  may  be  expressed  by 

- t x (t  x Vc)  £ Vc  - t t . Vc  - n n • Vc  . (1) 

It  is  this  component  of  Vc  which  causes  the  differential  progression  of  the 
wave  front  and  the  associated  bending  of  the  ray. 

Let  the  parameter  s be  a linear  measure  in  the  progression  of  the  ray. 
In  time  Increment  6t  the  length  Increment  traversed  is 

6s  * c 6t  . (*-) 

[n  terms  of  the  geometry  of  Pig.  I the  bending  of  the  ray  may  be  expressed 
by 


^ (3) 

6t  dn 

The  corresponding  governing  equation,  expressed  in  vector  form,  is 


— --  * t x ( t x Vc)  . (4) 

Dt 

The  differential  ope  ator,  D/Dt,  Implies  differentiation  in  time  moving  with 
the  ray.  Equation  (1)  readily  follows  from  Eq.  (3)  in  that  6_t  = - n 6a. 

The  change  in  a unit  vector  is  always  directed  normal  to  the  vector. 

It  the  speed  of  sound  is  relatively  constant  in  time  over  the  period 
of  propagation  of  the  ray  we  may  choose  to  integrate  along  path  length 
rather  than  in  terms  of  time.  In  this  case  Eq.  (4)  may  be  expressed  by 
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t x ( t x VK)  , 


whore 


K s In  c . 


2.  Vertical  and  Horizontal  Components 

Wo  proceed  by  relating  the  vector  equation,  Eq.  (5)  , to  the  natural 
reference  orientations  of  horizontal  and  vertical.  The  unit  vector  t , 
defined  as  before  as  the  local  direction  of  the  ray,  and  the  unit  vector  k , 
defined  as  the  direction  of  local  gravity,  together  generally  define  a 
vertical  plane.  This  plane  is  the  plane  of  Fig.  2. 


j x t 


Fig.  2 Vertical  and  Horizontal  References 


The  horizontal  component  of  t defines  the  orientation  of  the  unit 
vector  ^ . The  unit  vector  j completes  a right-handed  triple: 


k x i 
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And  j points  out  of  the  plane  of  Fig.  2,  toward  the  reader.  Whereas  the 
vector  k is  locally  defined  (i.e.  a function  of  position  only)  the  vectors 
j and  i depend  on  the  particular  ray  in  that  i is  parallel  to  the  horizontal 
component  of  t . The  angle  "Y  is  defined  as  shown  in  Fig.  2.  The 
parameter  Z denotes  the  depth  below  sea  level  of  the  location. 

Directional  change  in  the  unit  vector  t along  the  ray  path,  s, 
involves  only  two  components,  and  these  lie  in  the  plane  of  the  wave  front. 
The  plane  of  the  wave  front  may  be  defined  by  the  vector  pair  j and  j x t . 

The  vectors  _t  and  j x t may  be  expressed  by 
t = i cos  "Y  + k sin  y , 

j x t = - k cos  y + i sin  T . (9) 

We  proceed  by  developing  the  left-hand  side  of  Eg.  (5): 


KE  n 

— — = — ( i cos  y + k sin  V) 

D s Ds  — ~ 

ny  Di  Dk 

= (-  i sin  y + k cos  >)  ~ + cos  V — — + sin  V — — 

- - Ds  Ds  Ds 

Dy  , D i D i D k 

= - 0*  t)  — + cos  y (n  • + k k • + siny 

ny  D i Dk 

= - ( j x t)  — + cosy  i j,  • + (-  cosy  k i • + sin^)  Ds”.  10^ 

We  approximate  the  earth-curvature  factor  by 


cos  y 


Dk 

d7 


1 


R - z 


(ID 


where  R is  the  local  radius  of  curvature  of  the  earth.  Substitution  of 
Eq.  (11)  in  Eq.  (10)  yields: 


D t 


D7 


D» 


-os  y 


Ds 


r = - ( i x t ) ^ + cos  V u . — - ( J x t ) — _ z . 


02) 


Making  use  of  Eq.  (12)  we  now  take  the  ( j x t)  component  of 


Eq.  (5): 


D> 

Ds 


cos  y 


R-  z 


( j x t ) • t x ( t x VK) 


cos  y 


R - z 


+ ( j x t)  • VK 


cos  y 


R - z 


cos>  f + sin  y i . VK 
dz 


(13) 


And  the  j component: 


cos  y j 


D i 

dT 


. t x ( t x VK)  = - j • VK 


(1-0 


Equations  (13)  and  (14)  are  the  desired  pair  of  components  of  the  governing 
equation,  Eq.  (5). 

Provided  the  ray  does  not  become  vertical,  we  may  also  use  a linear 
measure,  x,  along  t , as  the  Independent  argument  for  the  integration: 


Dx  = Ds  cos  y . 


(15) 


We  may  also  introduce  the  linear  measure  y along  j . With  these  definitions, 
Eqs.  (13)  and  (14)  may  be  expressed  by 
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J A 


? K v 3K 

+ tan  > — 

?z  ox 


2 v SK 
- sec  ; ~ 
3y 


3.  Horizontal  Coordinates 

Let  c and  a'  be  a suitable  pair  oi  orthogonal  coordinates  for  the 
horizontal  expanse.  The  ascendents  define  the  corresponding  unit  vectors 
I and  I : 


E I , 


W I . 


We  choose  that  1 and  I form  a left-handed  system  with  k : 

I x I . - k . (20 

At  an  arbitrary  location  the  unit  vectors  i and  j , defined  eailiei  , 
form  angle  ft  with  l as  shown  in  1 ig.  3. 


rig.  3 Definition  of  the  Angle  fi 


We  proceed  as  follows: 


D i 


D x 


D 


- = j • — (cos  0^  + sin  0 J ) 


j • (-  sin  0 I + cos  0 I)  ^ + cos  0 + sin0 


123 


Di 


D0  + 

DX 


(cos  0 j - sin  0 j x k) 


Dx 

DI 

dT 


.DJ 

Dx 


P£ 

Dx 


DI 

Dx 


D£ 


Dx 


- J - ( i • v I ) . 


(21) 


Substitution  in  Eq.  (17)  yields 


m. 

D x 


Q + sec“  y ?K/?<y  , 


(22) 


where 


Q - J • (i  • VI)  (23) 

is  due  to  the  curvature  of  the  horizontal  coordinates.  Equation  (16)  is  the 
other  of  our  pair  of  governing  equations. 

•1.  Longitude  and  Latitude  Coordinates 

Longitude,  8 , and  latitude,  o , form  a pair  of  orthogonal  coordin- 
ates'*. The  ascendents  define  the  unit  vectors: 

Longitude  is  measured  east  of  Greenwich  only. 


V 0 * 

I / (R-z)  cos  O 

(2-1) 

Vo  3 

I / (R-z)  . 

(25) 

These  expressions  are  singular  at  the  poles.  The  curvature  term,  Lq.  (23), 
develops  as  follows: 

Q = - ) • | cos  f)  I + sin  fj  I | • VI 

» - J cos  3 1 • V I + sin  I • ^ I ^ 

- - j • cos  (J  J - cos  p tano/R-z.  (26) 

The  governing  equations  for  the  ray  trace  may  be  written: 


12?' 

*1 

-I 

Px 

R-Z  3z 

P£ 

cos  3 tan  O 

Dx 

R-z 

VK 


(3  7) 


(28) 


The  directional  derivatives  may  be  expressed  by 


I • VK  = 


J • VK  - 


1 


_ SK 
(R-z)  cos  o 38 

1 


R-z 


3R 

30 


(29) 


(30) 


Also  required,  as  increments  of  integration,  are 

6 8 = 7P  f*1--  6x  , 
(R-z)  cos  O 


(31) 


6 O = J 

sin  $ 

izr  6x  ■ 

(32) 

where  fix  is  true  horizontal  trace-length  increment. 

5.  Cartesian  Coordinates 

in  a Polar  Stereographic  Projection 

The  coordinates 

X 2 

L cos  0 

(33) 

y 5 

L sin  0 , 

(34) 

where 

L s 

(R-z)  M cos  p 

(35) 

M 2 

1 + sin  <pQ 

(36) 

1 + sin  c 

map  into  Cartesian  coordinates  in  a polar  stereographic  projection,  true  at 

latitude  o • The  coordinate  origin  is  at  the  pole, 
o 

The  ascendents  define  the  unit  vectors: 


vx 

s M l , 

(37) 

VY 

£MI. 

(38) 

The  curvature  term,  Eq.  (23),  transforms  as  follows: 
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I • ^ cos  0 I . V J + sin  0 I • V j j. 


cos  /3  sin  0 + sin  0 cos  0 


[ - sin<p 
z)  cos  o 


This  can  be  shown  by  expressing  the  unit  vectors  in  terms  of  the  polar 
vectors,  Eqs.  (24)  and  (25).  Elimination  of  0 and  <p  from  Eq.  (39),  using 
Eqs . (33)  through  (36),  leads  to 


sin  8X  - cos  0Y 
^ ~ (R-z)  “ (1  + sin  o ) 


The  governing  equations  for  the  ray  trace  may  be  written: 

& = - — - — + tan  y {cos  0 Mw  + sin  0 M ~ l 
Dx  R-z  dz  3X  °Y  ) 


D0  _ sin  0 X - cos  5 Y 
Dx  (R-z) “ (1  + sinoQ) 


+ sec2  y j sin  0 M ~ - cos  0 M 


6.  Integration  Formulae 

We  choose  to  carry  on  with  the  coordinates  introduced  in  Section  5: 
Cartesian  coordinates  in  a polar  stereographic  projection. 

The  governing  equations,  Eqs.  (41)  and  (42),  may  be  developed  for 
numerical  integration  by  introducing  finite  increments: 


- -{ih  + ff]6x  + tanr(f*  6X 


X 6Y  - Y6X 
(R-z)*"  (1  + sin  O ) 


+ sec  y 


3K  * V 

+ SY  6Y 


«-!y  6x 


To  these  we  may  add: 


6 X = 

M 

cos 

/3  6x 

(4  5) 

6 Y = 

M 

sin 

jS6x 

(46) 

6 Z = 

tan 

>6x 

(4  7) 

For  integrations  in  which  the  trace  is  permitted  to  become  vertical 
we  must  reintroduce 

6x  = cos  y 6s  , (-18) 

and  use  6s,  the  linear  progression  of  the  ray,  as  the  argument  of  inte- 
gration. 

For  integration  with  time  as  the  argument  we  reintroduce 

6s  = c 6t  . (-49) 


In  order  to  succintly  handle  surface  and  bottom  reflections  the 
argument  may  be  geared  to  depth  increment  6z,  unless  an  upper  bound 
on  6x  is  exceeded  in  which  case  6z  must  be  reduced. 

At  the  surface,  and  at  the  bottom,  reflections  are  performed  by 
reversing  the  sign  of  F.  At  the  bottom  this  involves  the  approximation 
that  the  bottom  is  flat. 


-11- 


Part  B:  A Scheme  for  the  Numerical  Integration 


We  herewith  develop  a proposed  integration  scheme  based  on  the 
geometric  visualization  of  refraction  and  the  adjustments  for  earth  and 
coordinate  curvatures.  Equations  (41)  and  (42)  express  the  ray  propagation 
as  differential  equations  in  terms  of  Cartesian  coordinates  in  a polar 
stereographic  projection.  In  order  to  visualize  the  geometry  we  rewrite 
these  equations  as  follows: 

curvature  term  refraction  term 


67  = — + j x t • VK  6s  (50) 

K-z  ~ — 

„ sin  /3X  - cos  BY  , 

6/3  = ^ K 6x  + sec  7 j . VK  6s  . (51) 

(R-z)  (1  + sin  cp  ) 

o 

The  vectors  j x t and  j may  be  recollected  from  Fig.  2:  The  vector  j x t 
lies  in  the  vertical  plane  and  is  normal  to  the  ray.  The  vector  j lies  in  the 
horizontal  plane  and  is  normal  to  the  ray. 

A segment  of  the  ray,  in  which  segment  the  refraction  terms  are 
taken  as  constant,  describes  a circular  arc  in  space.  We  propose  to 
integrate  the  ray  along  chords  for  such  arc  segments.  Representative  values 
for  the  refraction  terms  are  each  numerically  evaluated  for  the  segment. 

Along  a chord  the  curvature  terms  may  also  be  considered  steady. 

This  is  apparent  for  the  term  -6x/R-z  if  we  neglect  the  minor  local 
variation  of  R-z.  For  the  other  curvature  term  we  introduce  the  argument 
along  the  chord  increment  6s.  That  is, 


X = X + cos  0 7 
o 


(52) 


Y = Y + sin  p 7 . 
o ^ 


(53) 
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Substitution  in  the  coordinate  curvature  term  yield 


sin  gX  - cos  <3Y 
(R-z)2  (1  + sin  <p^) 


sin  /3  X o - cos  /3Y0 
(R-z)2  (1  + sinoo) 


(54) 


and  shows  that  this  term  may  also  be  considered  steady  along  a chord. 

We  assume  that  K,  the  natural  logarithm  of  the  speed  of  sound,  is 
given  at  an  array  of  grid  points  in  the  polar  stereographic  projection,  for 
each  of  a number  of  depth  levels. 

We  shall  describe  the  integration  in  terms  of  the  incremental  step 
from  position  Xr  , Yt  , z T , at  which  point  the  ray  is  directed  at  angles 

y ft  .to  the  next  position  where  these  values  are  indicated  by  sub- 

T T 

scripts  r + 1 . 

Let  y ^ . B ^ be  the  R-th  estimate  of  the  angles  with  which 
T+l  T+l 

the  ray  is  directed  at  position  r+  1.  The  first  guess  of  these  angles  may  be 
approximated  by 


y (o) 

T + 1 


= y_ 


A 


(o) 


T+l 


- A. 


(55) 


(56) 


Linear  extrapolation  from  T-  1 and  t values  may  yield  a better  first  guess. 
The  chord  is  directed  from  position  t with  the  angles 


(R) 

\ 

2 

(yT  ♦ 

y N) 

(57) 

(R) 

l 

= 2 

(sr  + 

ft  (R) 

^T+l  ) . 

(58) 
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Choosing  a chord  length,  6s  , the  R-th  estimate  ot  tin*  iny's  advance  ts 
given  hy 


6 X(K) 

•*  M 

JR)  , v (R)  . 

COS  p COS  y 6b  , 

(59) 

6Y^ 

* M 

. JR)  „ V(R)  . 

sin  p cos  > 6s  , 

(tiO) 

6z(,<> 

, y (R)  c 

■ sltr  7 6s  . 

(t.l) 

where  M is  an  appropriately  centered  value  ol  the  mapping  tactoi  , I'g.  (>U<)  . 
We  also  require 


6x 


(K) 


COS  > 


(K) 


6 s 


(lC) 


We  must  now  calculate  the  bending  experienced  by  the  lay,  ovei  this 
increment,  as  expressed  by  I'.qs.  (SO)  and  (SI). 

The  curvature  terms  are  evaluated  by 


R-z 


(t»d) 


o(R)  JR)  - (K) 

sin  8 X cos  8 Y 6 \ 

T T 


(R-z)' 


(1  + sin  $ ) 
o 


where  an  appropriately  centered  value  of  R-z  ts  list'd. 


(t>D 


Rot  evaluating  the  amount  ot  letractlon  wo  lequlie  lopiosontative 
values  of  the  refraction  terms.  We  propose  that  this  be  done  by  llnite- 
dlfference  evaluation  of  the  requiit'ti  directional  derivatives  ot  K.  l'he 
space  Increment  for  evaluating  those  should  be  commensurate  with  the 
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1 


length  of  the  chord;  In  fact  we  propose  that  it  be  made  equal  to  6s  to  give 
the  desired  representativeness. 


; 


B 


vertical  plane  horizontal  plane 

rig.  -1  Positions  A,  B,  R and  L foi  Evaluating  the  Directional  Derivatives 
of  K 


if 


The  locations  at  which  K should  be  interpolated  are  indicated  in 
rig.  -!.  In  tin'  vertical  plane  the  positions  are  A,  above,  and  B,  below; 
their  separation  is  6s  , forming  a cross  with  the  chord.  In  the  horizontal 
plane  the  positions  are  R,  right,  and  L,  left;  their  separation  is  6s, 
forming  a cross  with  the  chord.  The  R-th  estimates  ot  these  positions  are: 


X 

o 


J 

2 


(cos  } 


(R) 


v(R)r  /)  ,,  s 

sin  y ) cos  fi  M 6s 


I 


xJR)  - X 

B o 

(cos  , (R> 

- sin  x(R)) 

«(R) 

COS  P 

M 6s 

v « . V 

B o 

1 

+ — 

2 

(cos  V (R) 

- sinV<R)> 

A (R) 

sin  0 

M 6s 

.»  - « 

B o 

1 

+ 2 

(sin  V (R) 

* cosr(R)) 

6s 

x(R)  = X 

1 

2 

(cos  r(R) 

cosd(R)  ♦ 

sin/8) 

M 6s 

Y(«  = Y 

R o 

1 

+ 2 

(cos  r(li> 

sin  S(R>  - 

cos  „(R)) 

M 6s 

(R) 

1 

>-  — 

Sin  V <« 

6s 

ZR  = Zo 

2 

x<R)  * X 

L o 

♦ i- 
2 

(cos  r(R) 

a(R) 

cos  0 

sin  d(R)) 

M 6s 

V,(R)  - Y 

L o 

1 

+ — 

2 

(cos  r(R) 

sin  0 (R)  * cos  A 

M 6s 

(R) 

z . = z 

L o 

1 

+ — 

2 

sin  > <W 

6s  . 

(65) 

Let  the  values  of  K interpolated  at  these  locations  be  denoted  by 

..  (R) 
ka 

, K 

(R) 

B 

K,f»  and 

w (R) 

kl 

• 

(66) 

refraction  terms  are 

accoulinyly  given  by 

( ) x t 

• V K 

6s)  (R> 

..  (R) 
ka 

k(r) 

kb 

(67) 

see  VlR)  ( , 

. « Ss)  M - 

sec  >(R)  (K 

(R)  ..  (R) 

R " \ 

) 

(68) 
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The  R-th  estimate  of  the  bending  of  the  ray  in  the  coordinate  system  is 
given  by 


R-z 


+ 


(69) 


sin  |9(R)  X - cos  /3(R)  Y ,d. 

__2 _ 2 6x(R) 

(R-z)“  (1  + sin  <p^) 


sec 


- K ‘iy> 


(VO) 


These  give  the  R+l  estimate  of  the  angles  with  which  the  ray  is  directed 
at  position  r+l: 


I 


y (R+l) 
T+l 


=3 


+ <>y 


(71) 


« (R+ 1) 

^T+ 


i 


0. 


(72) 


With  these  new  estimates  the  R-l  cycle  is  performed  beginning  again  with 
Eq.  (57).  The  cycle  is  repeated  until  convergence  is  attained  (i.e.  within 
some  finite  tolerance).  Upon  such  convergence  the  desired  chord  iias  been 
found,  the  incremental  step  is  completed,  and  the  ray  moves  to  position 
T+l. 

It  should  be  noted  that  in  evaluating  the  coordinate  curvature  term, 
Eq.  (64),  the  origin  for  X and  Y is  at  the  pole.  Also,  the  angle  is  the 
angle  made  by  the  ray  with  the  X axis;  that  is, 

cos  0 ■ • .1  . (73) 


L 
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In  the  suggested  scheme , the  chord  length,  6s  , may  be  chosen 
for  each  step.  We  teel  that  the  optimum  choice  tor  6s  would  be  such 
that  convergence  ts  attained  at  R .1.  However,  it  may  be  worth  experi- 
menting with  6s  to  see  how  the  convergence  is  affected. 

Reflection  ot  the  ray,  at  the  sea  suitace  and  at  the  bottom  also 
influences  the  choice  ot  6s.  l.et  the  distance  to  the  boundary  being 
approached  by  the  ray  be 


6/.* 


(74) 


where  Z is  zero , or  the  depth,  h.  It 


then 


6 s 


(K»  1) 


sin  V 


(K) 


(7!>) 


should  be  used  as  chord  length  in  the  next  cycle,  so  that  the  step  takes  the 
lay  exactly  to  the  boundary . At  this  juncture,  Y ts  replaced  by  -Y , 
assuming  the  bottom  Is  horizontal . 

Another  contingency  must  be  provided  tor.  In  approaching  the 
upper  ot  lower  boundary,  one  ot  the  other  ot  the  A,  above,  or  It,  below, 
locations  may  tall  beyond  the  boundary.  Provision  should  be  made  so  that 
the  interpolation  of  K will  yield  a suitable  extrapolated  value  in  such 
cases . 
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Flow  Diagram  in  Terms  of  Operations 


Position: 
Angles  of  Ray: 


X Y Z 

T T T 

VT 


Set  estimates:  Eqs.  (55)  and  (56) 


Compute  chord  angles:  Eqs.  (57)  and  (58) 


Choose  6s 

Is  6z  too  large?  No;  proceed  with  6s. 

Yes;  replace  6s  according  to  Eq.  (75). 

Compute  curvature  terms:  Eqs.  (63)  and  (64) 

Compute  locations  A,  B,  R,  L:  Eq.  (65) 

Interpolate  K at  A,  B,  R and  L 

Compute  refraction:  Eqs.  (67)  and  (68) 

Compute  total  bending:  Eqs.  (69)  and  (70) 

Compute  y and  8 : Eqs.  (71)  and  (72) 

Test  angles  for  convergence: 

No;  recycle 

Yes;  advance  step  (including  any  reflection)  to  T+  1 
position  and  angles 


The  program  requires  an  interpolation  subroutine  for  interpolating  to 
specific  locations  in  the  three-dimensional  distribution  of  K = In  c,  where 
c is  the  speed  of  sound.  Computationally,  the  program  requires  zero-order 
continuity  of  the  interpolation;  hence  even  linear  interpolation  between  grid- 
point  values  could  be  used.  However  in  the  interests  of  better  accuracy 
higher-order  interpolation  is  preferable. 

The  three-dimensional  distribution  of  K may  be  expressed  by  the 
two-dimensional  distributions  for  a number  of  discrete  depth  levels.  These 
depth  levels  should  be  concentrated  where  the  physically  and  operationally 
significant  variability  of  K is  greatest.  The  set  of  fields  is  represented  in 
terms  of  a polar-stercographic  grid-point  array.  It  is  quite  feasible  that 
the  three-dimensional  distribution  of  K,  for  an  ocean  region  of  interest, 
could  be  stored  in  the  core  memory  of  the  CDC  6400  computer. 

The  interpolation,  for  location  XA  < YA  < ZA  for  example,  may  be 
performed  by  interpolating  to  at  the  four  levels  which  include  the 

depth  z^  in  the  middle  layer.  Cubic  interpolation  in  z would  then  com- 
plete the  interpolation.  It  should  be  recalled  that  z^  may  lie  above  the 
sea  surface;  in  this  case  linear  extrapolation  of  K from  the  surface  and 
first  sub-surface  depth  may  be  adequate. 
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